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ABSTRACT 


This thesis addresses the problem of identifying the 
dynamics of a linear system in the frequency domain. An 
alogrithm operating on the Fast Fourier Transform (FFT) of 
blocks of signals is developed and its performance evaluated 
through computer simulations. Several properties are tested, 
in particular, its convergence and its capabilities of 


identifying the frequency response of the unknown system. 
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I. INTRODUCTION 

The problem of identifying the dynamics of a physical 
system from its input-output behavior is a fundamental issue 
in modeling and control. The most general approach is based 
on the assumption of a mathematical model, often linear and 
time invariant, and its parameters are determined by 
optimization of an error criterion. 

A vide class of identification algorithms exists in the 
literature. Several books have been vritten on this matter, 
see for example (Ref. 1, 21. Of particular interest are the 
recursive algorithms, vhere the estimated parameters are 
updated on-line based on nev input-output data collected at 
each sampling instant. In most of the approaches on recursive 
identification the parameter estimates are updated for every 
new set of data points; as a consequence, the estimate at any 
time is updated up to the current data. 

A different approach is based on Block Processing (BP), 
where the parameters are recursively updated on the basis of 
the data within nonoverlapping e EE time. Although 
this introduces a delay between the time we measure the input- 
output data and the time we use it in the estimation, it has 
the advantage that time averaging within the time interval 


reduces the effect of disturbances. 


Block Processing techniques have been investigated by 
several authors. Of particular interest for this thesis is 
the approach developed by Mansour [Ref. 3] in the adaptive 
filtering contest, where the parameters are updated on the 
basis of the DFT (Discrete Fourier Transform) of the data 
wlthin each time interval. 

Im this thesis an adaptive identification algorithm 
operating on the DFT (Discrete Fourier Transform) of blocks of 
input and output data is introduced. In particular, the 
algorithm identifies the frequency response of the system's 
linear model by applying several recursive estimation 
techniques and adapting them to the frequency domain approach. 
We will investigate estimation methods based on the Projection 
Algorithm and Recursive Least Squares [Ref. 4] and the results 
will be compared. Even though Recursive Least Squares is a 
complex and time consuming method, it is preferable when 
accuracy and speed of convergence are needed. 

This thesis is organized as follows. Chapter II 
introduces the concept of Block Processing and develops two 
algorithms for recursive identification in the frequency 
domain. Proof of convergence is also part of this chapter. 
Simulation studies are introduced in Chapter III showing the 
effectiveness of the algorithms and the conclusion is in 


Chapter IV. 


II. IDENTIFICATION OF LINEAR SYSTEMS 


A. STATEMENT OF THE PROBLEM 

Consider an LTI (Linear Time Invarlant) system with input, 
output signals given by u(t) and y(t) respectively. The 
problem we address is the estimation of the parameters of a 
linear model in order to predict future values of the output 
y given past measurements of the input and output signals. 
The general block diagram is given in figure 2.1. The block 


211 denotes time delay by one clock pulse. 





Figure 2.1 Block diagram of the system 


This figure shows the general scheme where y(t) is the 
prediction of y(t) based on the past input, output, and the 
estimated model. The purpose is to determine a model for the 


plant so that the error e(t) is miniminized. In particular, 


the model is restricted to be an LTI system in either one of 
the two forms: 


ARMA (Auto Regressive Moving Average) 


y(t) =8,y(t-1)+...+28 PEA m Da 


where a,, b, are constants; or 
CONVOLUTION form, where the output input signals are related 


by the convolution sun. 


y (ah IN (2.2) 


y(t) Y. h(x)u(t-x) (250 


T=-00 


The problem of adaptive frequency response identification 
mainly has been concentrated on the approximation of infinite 
impulse response systems by finite impulse response models. 
In this thesis, we will consider the identification of stable 
plants only, for which the impulse response h(t) decays to 
zero as time increases. For this reason, the infinite impulse 
response h(t) can be approximated by a finite sequence with 
length N. In this way, the moving average (MA) model of the 


following form is obtained. 


y(t)=h(O)u(t) +h(1)u(t-1) +-+h(N-1)u(t-(N-1)) (2.4) 


where h(t), for 0 < t < N-1, is the truncated impulse 
response, We choose the block length N to be a power of 2 in 
order to apply FFT (Fast Fourier Transform) techniques. 

The approach we consider is based on a frequency response 
framework. In particular, the estimated model is computed 
from the Fourier Transform of blocks of data collected from 
measurements of the input and output. 

There are several reasons which motivate this approach. 
In many cases of interest, a model which is accurate on a 
desired frequency range only is identified for the plant. For 
example, 1f we want to model the low frequency spectrum of the 
system with the frequency response approach, we select weights 
in order to select frequencies of interest. 

In the next section, a recursive algorithm for the 


identificaton of the frequency response is introduced. 


B. RECURSIVE FREQUENCY DOMAIN APPROACH 

Block Processing (BP) techniques for on line 
identification of linear models have been investigated by 
several authors [Ref. 3]. With this approach the estimated 
model parameters are sequentially updated on the basis of 
blocks of data, rather than at each data point. In this 


section, we investigate a BP adaptive algorithm which operates 


on the Fast Fourier Transform (FFT) of blocks of input and 


output data. 


blocks of data blocks of data 


——“— 


4N-1 ON-1 


mu. 


4N-1 ON-1 





Figure 2.2 Blocks of Input and Output Data 


There are some advantages in dealing with data by blocks 
rather than individually. Block processing averages the data 
within the block, resulting in better noise rejecton. In 
order to illustrate the BP approach in the ARMA model case, 
consider the model in (2.1) and the data y,(t), u,(t) within 


the k-th time block, i.e. 


yx(t) zy (KN+t) 


u,(t) =u(kN*+t) (2.5) 


for t = О, 1, ... , N-1. If the block size N is larger than 


the system order n, we can write (2.1) as 


pM 7007 Tsay (t-2)+b u, (t-1)+b.u,(t-2) (216) 


Mona Second order case with 2 < t < N-1. It is easy to see 


that (2.6) can be written in convolution form as 


yx(t) za(t) Gy, (t) +b(t) Bu, ( t) n<t<N-1 (2.7) 


with a and b sequences in R^ 


a=[0 a, a, ase a, Ди 0) 


b=[0 b, b, = b, 0 — 0] (2.8) 


and @ denoting circular convolution. 

In the approach, the well known relationship between the 
DFT and the circular convolution in the frequency domain is 
exploited. By defining { Y,(1) , 1=0,1, ...,N-1 ) and ( U,(1) 
, 1=0,1, ...,N-1 ) as the DFT of the output and input data 


y, (€) ,u,(t) in the k-th block, we can define Y (1) as 


AI AMO A 1-0,1,-,N-1. (2.9) 


Clearly from (2.7) and (2.9), we can write the important 


relationship 


yy (£) =P, (€) nstsN-ı (2.10) 


with Y ,(t) =IDFT[ Ç ,(1)], and IDFT denoting inverse DETE 
A linear estimation algorithm for A(1) and B(1) in (2.9) 


can be determined by defining an error signal 


d t) -$ t<N- 
where $ ,(t) =IDFT{ A ,(1)X (1) +B (1)U (1), 1=0,1,...,N-1 ] and 


А,, B , are estimates available at the k-th block prong. 


above definitions, we can write the recursive estimates for A, 


B by the algorithm. Let 


Ә, (1) =|А,(1),В,(Ј)17 


ANA Gl s (2.12) 
Гог 1 = О, 1, ..., N-1, k cx. 07" SU ə Ve Le 
E, (1) =DFT[é, (t)] (2.13) 
with é, as in (2.11). Then the recursion is such that 


5:75: 
6... (2) “6. (2) 4 — -- (2.14) 
¿max 
DPD) 


For this algorithm we can show the following: 
EEG. 1) S ©41) |“ < | 6,(1) - (1) 1“ (2.15) 
for all 0 1, 2, ..., N-1, and for all k: 


(b) If the input and output data u,(t), y (t) are bounded for 
all k and t, then 


T e,(t)=o for all OstsN-1 (2.16) 


In this result we see that the parameter error decreases 
with time (eqn. (2.14)) and the error between the model and 
c plant e (t) tends to zero, 1.e. the model output tends to 
follow the plant output. 


Proof: Let us define the parameter error 


6,(1) -6,(1) -0,(1) (2.17) 


Then from (2.14) and (2.17) we can write the recursion 


6,, (1) 26,01) -u, Ga (1) Ë. (1) ) (2.18) 


where 


1 


k 1 (2)) 2+|Y, (1) |) (219) 


By taking the magnitude square of both sides of (2.18), we 


obtain 


(Gran (1) 261,118, (1) 
-[6'7, (21) -u E" 7, (1) Xe (2)] Ic CD nox (1) EL (2) 1 


=|6,, (1) |? Bg? (1) Xe (1) 6, (1) - 2077 (1) Xe (1) ÈD) 


+ppe, Xp (1) Xx (1) B, (1) (2220) 


By taking the sum of all frequency components in both sides of 


(2.20), we obtain the expression 


N-1 


E N-1 
Nr (DP 3: B. CD O E UA 0 (1) 
1=0 1=0 1-0 


N-1 


Y BAD AR ADE I) 
1=0 


N-1 


> WERE ER (21209) 
m0 


Applying Parseval's theorem, ve can relate data in the 


frequency domain to the corresponding time domain as 


10 


N-1 
E AO MB) 6 (tie k ©) 


1-0 t=0 
N-1 N-1 
De Oea Eeo 
t=0 t=2 


where we define 


Ey (E) =IDFTI[X: (1)6,(1)] 


Equation (2.23) comes from the fact that 


ë, (t) =0 (SÈ AZI 
6560706 6) n<t<N-1 
and 
e, (€) =e,(€) N] 


Still using Parseval's theorem, we can write 


N- 


Y p, Copt pu 


1-0 


in 


os 


(2 


(2. 


e. 


(22 


(2. 


(2. 


2 


23) 


24) 


25) 


26) 


27) 


28) 


After some manipulations and the fact that j | X (1) j° < 1 fen 
all k and 1 = 0, 1, ..., N-1, we can bound the parameter error 


at the end of the (k+1)th block as 


Pro: 








° < Bd 2 ëd +u ëd 2 


and therefore 


Px. s 4^ 71484 (2.30) 








Finally, since || e, | 2 0 for all k and is a nonincreasing 
sequence as in (2.30), the increment DN e li must tend to zero 
as k tends to infinity. Therefore 


lim I”A” 
k 


Soot -0 (210989 
1*MAXIY, G2) 0, CD P] 


which proves the result. 


The significance of this result is that we can estimate 
the parameters of the given system (in terms of Á (1), B (1) 
based on the frequency content of each block of data. 
Furthermore, the estimation is recursive, and the convergence 
of the prediction error e(t) to zero is guaranteed at least in 
the ideal case. It will be shown in the simulations that even 
in the presence of measurement noise, the convergence of the 


error to small values is still satisfactory. 


rz 


In the implementation of the estimation algorithm 
introduced above, we need to compute the DFT of the "windowed" 
nem term e, (t) in (2.25) and (2.26). We can see that in the 
implementation two operations of IDFT are required to compute 
Y, (t). Using DFT again, E,(1) is obtained, which is used in 
the recursion (2.14). 

An algorithm, which does not require this sequence of IDFT 
and DFT in its implemention, is introduced next. The new 
algorithm also offers the advantage of being able to use 
different recursive estimation techniques (such as the 
Recursive Least Squares) which exhibits faster convergence. 
However, this is obtained at the expense of added complexity. 

in order to introduce the technique, let us recall the 


signal 


¥,(t) =a(t) Sy, (tc) +b(t) Wu, (Cc) (2.32) 


with a, b being the plant parameters, and the equality 


Yg (t) =; (t) nstsN-1 (2.33) 


Going to vector notation, let us write (2.33) in the following 


form 


43 


0 Ӱ, (0) 


= ак d [a [LC] zi p (n) aa 
- Do) Iz] .. 1 
y,(N-2) Ee 
bal ig 





where [0] and [I] denote blocks of zeros and the identity 


matrix respectively. Definition of the Fourier matrix F as 
7— 
o (2.35) 


allovs the computation of the DFT of a sequence as a matrix 


operation. In this way, we obtain from (2.34) 


Ру,=Еһӱ, (2236) 


and therefore, after simple manipulations, the following 


equation, 


Fy z (FhF“) (FY,) (2.37) 
It is easy to see that the vectors Fy, and F y y, are the arrays 


of DFT coefficients of the respective sequences y, and Y. 


which yields 


Recall from (2.32) that 
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Ý, (1) 2A, (1) V, (1) +6, (1) U, (1) (2.39) 


for 1-0, 1, 2, ..., N-1. Let us define YU as the following 


expression 
mio) 0 = 0 AU (O O 
0 a ace О О D = 0 
YU= 5 k (2.40) 
0 e : : : 0 ` 


and combine (2.38) and (2.39) to yield 
A, (0) 
A, (N-1) 
Y YU s (2.41) 


B, (0) 


B, (N-1) 


The recursive estimation algorithm is derived from (2.41). In 
particular, we know that (2.41) can be broken into a sequence 


of linear equations of the form 


Y, (1) =0;(1)0 10575: (2.42) 


where 


15 


8=[A,(0) -- A,(N-1) B,(0) + B,(N-1)]7 (2.43) 


is the vector of unknown parameters, and 
$r(1)-[H(1,0) — H(1,N-1)] (YU) 
-[H(1,0) Y, (ODE IA 


¿H(1,0)U,(0) — H(1,N-1) U,(N-1)] 
(2.44) 


In (2.44) the coefficient H(i,j) are constant (in the sense 
that they do not depend on the block k) and are determined by 
the window matrix h only. 

From (2.42) we can use any recursive algorithm to estimate 
0. In particular, a Recursive Least Squares (Ref. 2] 
algorithm applied to complex data is going to yield the 
following recursion 


P®, (1) (Y, (1) -0%(1)0]) 


011-01. 
PES 1+Ô7(1) PLO, (1) 


(2.45) 
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P,o,.(1)Öx(1) PL 


Pi"-pl- - : (2.46) 
1*0,(1) PE ,(1) 
57:52 = 0, 1, 2, ..., N-1, with initial conditions at the 
beginning of each time block given by 
62-87. ; Pk=o2I (2.47) 


^ an arbitrary constant 


with I the identity matrix and o 
parameter. In general, the constant o is chosen to be a large 
value, in comparison with the order of magnitude of the 
parameters 6. 

Although far more complicated to implement, the algorithm 
(2.45), (2.46), and (2.47) has better convergence properties 
than the one shown previously. 


In the next section the two algorithms are applied to the 


identification of a linear time invariant system. 


1 7 


III. SIMULATION STUDIES 


A. INTRODUCTION 

The estimation techniques discussed in the previous 
chapter have been simulated, and used to estimate the 
frequency response of a systen. In addition, two computer 
programs written in MATLAB have been developed. The first one 
is the main program given in Appendix A which simulates the 
entire system in an iterative fashion. The main program 
estimates the value of the frequency response of a given 
system at a particular frequency and compares it with the 
corresponding value of the original system. 

The second program is a function type of subroutine which 
implements the Recursive Least Sguares algorithm [Ref. 2]. 

The proposed identification method has been tested using 
several different inputs in the program in order to simulate 
and analyze the effect of various excitaton conditions. In 
the next section, a first order recursive difference equation 


is used as an example. 


B. SIMPLE DERIVATION 
The system being used as an example is a first order 


system described by the difference equation. 


18 


y(t) =0.3y(t-1)+5.0u(t-1) (3010) 


in like manner (3.1) can also be written in terms of the 


impulse response, 


y (€) =h(0)u(0) +h(1)u(t-1)+.. 23 him) u(t-m) (5) 


mz -cə 


Since the system is stable, its impulse response decays to 


zero, and (3.2) can be approximated by a finite sum 


y(t) =h(0) u(t) +h(1)u(t-1)+...+h(N-1) u( t-N+1) (3.53) 


In the followig simulations N is chosen as N = 32 data points. 


The transfer function H(z) of the original system 








Y(z) 5 
Н(2) = = 3.4 
- U 5.50 ) 
yields its frequency response 
B 5 
Ба ыы 


Due to the use of FFT methods, we estimate the frequency 
response at discrete frequency points © = 2m71/N, 172 0, 1, 2, 


..., N-1. For this purpose, let us define 


129 


+00 T ELM 
H(1)*$ hine (3.6) 


n=0 





and its finite approximation 


j( 2.12) 


N-1 E 
#(l)=) h(n)e N (3.7) 
n=0 


The next step is to calculate the frequeccy response of the 


original system. 


C. STRUCTURE OF THE PROGRAM 

The algorithms introduced in the last chapter have been 
tested in the identification of a discrete time system. In 
particular we look at the problem of identifying the frequency 
response of the system by fitting a finite impulse response 
system to the input-output data. 

The software developed consists of two programs: a main 
program which produces the input output data used in 
identification, and a subroutine which implements the 
Recursive Least Squares identification shown in Chapter II. 

The effectiveness of the estimation is assessed on the 
basis of two criteria: the error between actual output of the 
system and predicted output, and the error between the 
frequency response of the system and the frequency response of 


its estimate. 
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Since we try to fit a finite impulse response model to a 
recursive system (with impulse response of infinite length), 
we cannot pretend to estimate the frequency response at all 
frequencies in an effective way. Therefore, we input test 
signals at finite number of frequencies and we look at the 
estimated spectrum at these frequencies only. 

As it will be seen in the sequel, the results are 
satisfactory even in the presence of added disturbances in the 
measurements. The recursive use of blocks of data has the 
effect of smoothing the effect of disturbances in the 


estimates. 


D. SIMULATION RESULTS 

In this thesis, ve try inputs vith different frequencies 
to test the identification algoritnhm. Also, we test the 
behavior when the input is the combination of several 
different frequencies to examine the effectiveness of the 
technique. In the next paragraph, we discuss some results 
obtained by running the program. 

1. Original System 

Figure 3.1 shows the magnitude and phase of frequency 

response of the original system. It shows the low pass nature 


of the system. 
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2. Single Frequency Input at Arbitrary Frequency without 
Noise 

Figure 3.2 and 3.3 are in the same group in which we 
tried a constant input with magnitude 10. As predicted, 
figure 3.2 portrays the error EPS between the output and 
predicted output which decreases and approaches zero. 
Furthermore, in figure 3.3, the frequency response of the 
estimated model at the input frequency is shown to be 
identical to the one of the original system. 

In a second set of experiments, we tried several 
sinusoidal inputs with different frequencies. First of all, 
we used an input signal with a digital frequency 7/2 and 
obtain the results shown in figures 3.4 and 3.5. And then, we 
tried another input with frequency 7/8 and get the results 
shown in figures 3.6 and 3.7. In like manner, the error 
approaches to zero quickly and the frequency response of the 
estimated model at the input frequency is the same as the one 
of the orignal system. 

3. Multiple Frequency Inputs without Noise 

In a third set of experiments, we tried the 
combination of four frequency inputs at arbitrary frequencies 
to obtain the results shown in figures 3.8, 3.9, and 3.10. 
The output in the time domain is shown in figure 3.8. In this 


case, the error between true and predicted outputs approaches 


22 


zero quickly and the frequency responses at the input 
frequencies are still identical. 
4. Single Frequency Input with Noise 

In this section we show results of similar experiments 
wlth random measurement noise. The noise considered is 
Gaussian white zero mean, l.e., lndependently and indentically 
distributed. Its standard deviation is set at about one third 
of the output magnitude. 

As before we first excite the system with a constant 
input, and we add measurement noise to the output of the 
system. The results are shown in figures 3.11, 3.12, and 
3.13. From figure 3.12, even though we find that the error 
does not decay to zero due to the presence of noise, the 
frequency response of the estimated model is sufficiently 
close to the one of the original system. 

In a second set of experiments, we tried sinusoidal 
inputs again with noisy measurements. There are four inputs 
with different digital frequencies (117/16,7/2, 7/4, and 7/16) 
used and the results are shown in figures (3.14-3.22) 
respectively. Then, by inspecting the figures, the error 
still does not go to zero due to the presence of noise, but 
the estimated frequency response is sufficiently close to the 


one of the original system. 
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5. Multiple Frequency Inputs with Noise 

In like manner, we input signals with multiple 
frequencies to test the alogrithm. In this case, we used four 
different groups of inputs. The first one is constant and 
sinusoidal (frequency 7/2) input. The 2nd one is a constant 
and sinusoidal (frequency 7/4) input. The 3rd and 4th are 
sinusoidal with two frequencies (7/2,7/4) and (7/4,7/8) 
respectively. The results are shown in figures 3.23 to 3.31. 
Then, inspection of these figures reveals that the error in 
each case is bounded and the frequency responses of the 
estimated model are sufficiently close to the one of the 


original system. 
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IV. CONCLUSION 

Im this thesis we introduced an approach to the 
identification of linear models based on the frequency domain 
formulation. The recursive nature and the use of FFT 
techniques make the approach attractive for on-line 
implementation. 

A particular aspect we addressed is the identification of 
the frequency response of the system by approximating the 
impulse response with a finite sequence. By use of a 
simulated example we have explored the convergence properties 
of the algorithm, and its robustness in the presence of 
measurement noise. 

Several issues remain to be investigated, mainly the 
advantages (if any) of this approach in comparison with more 
conventional time domain estimation techniques. Although this 
is subject of future research, as a conjecture we can say that 
the recursive frequency domain approach might be more robust 
than the conventional time domain in the cases when the input 
signal has energy concentrated around a few frequencies, such 
as the case of periodic excitation. In this case the 
processing gain of the DFT should give a better performance in 


the presence of measurement noise. 
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APPENDIX A : MAINPROGRAM 


THESIS MAIN PROGRAM File name : THESIS4.M 


Advisor : Roberto Cristi 
Student : Chao, Chi-Shun 
Date SAI 1091 


o9 o? o9 o9 o? o9 o9 o? o9 


oo 


Set some constant number 


elg 

clear 

tfinal = 80; 

de= 0.1, 

kmax = tfinal/dt; 
mana ('normal'):; 
rand_mag = 35; 

M = 20; 

N = 32; 

N2= N/2; 


The max. time to run 
Samplinç time 

The max. number of signals 
Set the type of noise 

Set the gain of noise 

The # to calculate the error 
The # of signals in a block 


oi oi oi o oi oi 40 


% Set the initial condition 
Ek (est 

yp(2)=0; 

paco; 

y (2)=0; 

u(l)=O; 

u(2)=0; 


Produce the input and coresponding output 
at different frequencies 
or n=3:kmax 


Fh o9 oe 


2000710” 
u(n)= 20*cos(pi*n*5/16) ; 
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oi oi oi SN AS QO QO QO 


30*cos(pi*n/4); 

40*cos(pi*n/8); 
50*cos(pi*n/16); 
60*cos(pi*n/32); 
10-30*cos(pi*n/4); 
40*cos(pi*n/8)+30*cos(pi*n/4); 
10+20*cos(pi*n/2)+30*cos(pi*n/4)+40*cos(pi*n/8); 
10+20*cos (pi*n/2)+30*cos(pi*n/8) : 

+40*cos (pi*n/8+50*cos (pi*n/16) ; 


OO oN o0 QO oO QO oe 
al: NE EE 
Sows Jeo 14 13 


ut Nut Yet Yet Yet ee Non Se” 
l Hu H H w ıı ı 


yp(n)=0.3*yp (n>1)+T5 ии hp; 
y(n)=yp(n)+trand_mag*rand; 


end 


% Produce the figures of input and output 
piot (ü); 

title( INPUT UTu2” 

xlabel('n'); ylabel('MAGNITUDE') ; 

grid: 

pause; 


cig 

subprot (2111F Pone 

title('YP(n) ORIGINAL OUTPUT WITHOUT NOISE '); 
xlabel('n'); yilabel ("MAGNITUDE 

grid, 


subplot (22), plot(y); 

title('Y(n) OUTPUT WITH NOISE u(n)=20cos(5*pi*n/16) 
yino p CIE; 

xlabel('n'); ylabel('MAGNITUDE!); 

grid; 

%meta thesis44 

pause; 
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% Produce the frequency response of the original system 
y using the "dbode" function 

num=(5]; 

den=[1,-0.3); 

w-linspace(0,pi,17); 

[mag,phase]=dbode(num,den,w) ; 


elg 

52 Plot(211),ploti(v,maq), 

title('FREQUENCY RESPONSE OF THE ORIGINAL SYSTEM H(Z)'); 
Zabel ('THETA'); 

ylabel('MAGNITUDE'); 

g id; 


subplot (212) plot(w, phase); 

title('FREQUENCY RESPONSE OF THE ORIGINAL SYSTEM H(Z)'); 
? cbel( THETA"), 

ylabel( " PHASE (DEGREE)!); 

amd; 

pause; 


% Create the matrix H 
v=[zeros(1,N2),ones(1,N2) ]; 
h=diag(v); 

F=fft(eye(N)); 
H=F*h*inv(F) ; 


% Set the zero matrix 
THETAK = zeros(N,1); 
PO=10*eve(N); 

P=PO; 
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% Use the FFT and RLS (Recursive Least Square) techniques 
% to estimate freguency response of the model. It is an 
% interactive and iterative loop 
for k = 2 M 
yko = y( (k-1)*N2+1 : k*N2 )'; 
uk — Uf. (KT21+N211 > 8A ee 
yk=[zeros(N2,1) ;yK0]; 


YK = fft( yk ); 
UK = fft( uk ); 
EPS (k)20.0; 


for 1S PN 
PHI=diag(H(i,:))*UK; 
[THETAK,P]-rls(THETAK,P,PHI,YK(i)); 
EPS(k)=EPS(k)+(YK(i)-conj (PHI') *THETAK) '* 
(YK(1) con] (PHI' |= FRE TAK 
end 
THETA(:,k) = THETAK ; 
end 
EPS = sqrt(EPS)/N; 
% Calculate the corresponding horizontal axis 
for k=0:N/2 


w1l(k+1)=2*pl1*k/N; 
end 


ee U. uu > E 
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% Produce the frequency response figures of estimated model 
clg 

20 .06(211)  plot(vl,abs(THETA(1:17,M))): 

555 -("EST. H(2) '); 

xlabel('THETA (65:777(1N/2) 0” 7 

ylabel('MAGNITUDE'); 

grid; 


1 5150t(212) plot(w1,angle(THETA(1:17,M))): 
title( "EST. H(2) '); 


xlabel ( " THETA ( PT=K/(N/2) )?); 
ylabel('PHASE (RADIAN)'); 

gend; 

pause; 


% Produce the error figures 

cig 

tlisplot (EPS): 

title( "EPS (when u(n)=20cos(5*pi*n/16) 
pa SEENEN E 

Blabel('K'); 

ylabel ( 'MAGNITUDE'); 

grid; 

%meta thesis44 

pause; 


Produce the frequency response of the estimated model 
in a continuous pattern. 

c TG 

57 .(V,mag,vl,abs(THETA(1:1/7,M))): 

title('H(Z) & EST. H(Z) '); 

xlabel('THETA ( PI*K/(N/2) )'); 

ylabel('MAGNITUDE'); 

grid; 

pause; 


A NO 
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% Store the value of the response of the estimated model 
% at the input frequency. 
THETA_F = zeros (17,1); 


%THETA F(1,1) = THETA(1,M); 
STHETA F(2,1) = THETA(2,M); 
$THETA F(3,1) = THETA(3,M); 
$THETA F(5,1) = THETA(5,M); 
$THETA F(9,1) = THETA(9,M); 
THETA F(6,1) = THETA(6,M): 


Produce the freguency response of the original system and 

the response of estimated model at the input freguency 

in the same figure. 

C LO 

plot(w, may, 

title('H(Z) & ESTIMATED H(Z) AT THE INPUT FREQUENCY ... 
(WITH NOISE)'); 

xlabel ( 'THETA ( PI*K/(N/2) )'); 

ylabel( "MAGNITUDE"), 

grid; 

hold on 

displot(wl,abs(THETA F)); 

Hora o me 

meta thesis44 


oo o9 o9 
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o9 o? o9 o? o? o? o? o? o9 o? o9 o9 o? o? o9 o? o9 o? o? o? o? o? o? o? o? 


ove oo oo oo 
ov oo oe oo 
ovo oo до — oo 
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